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Abstract 

We present a new spectral method for the Direct Numerical Simulation of Magnetohydro- 
dynamic turbulence at low Magnetic Reynolds number. The originality of our approach is 
that instead of using traditional bases of functions, it relies on the basis of eigenmodes of the 
dissipation operator, which represents viscous and Joule dissipation. We apply this idea to the 
simple case of a periodic domain in the three directions of space, with an homogeneous mag- 
netic field in the direction. The basis is then still as subset of the Fourier space, but ordered 
by growing linear decay rate |A| {i.e according to the least dissipative modes). We show that 
because the lines of constant energy tend to follow those of constant |A| in the Fourier space, the 
scaling for the the smallest scales |A™^^| in a forced flow can be expressed using this single pa- 
rameter, as a function of the Reynolds number as y^\X^^^\/{2nkf) :^ 0.5Re^^^ , where kf is the 
forcing wavelength, or as a function of the Grashof number C/, which gives a non-dimensional 
measure of the forcing, as |A™^^|^/^/(27r fef ) ^ 0.47G?-^°. Th is sc aling is also found con- 
sisten t with heuristic scalings derived by lAlemanv et al\ (|l979[ ) and IPotherat Alboussierd 
(|2003h for interaction parameter 5* > 1, and which we are able to numerically quantify as 
k^^^'/kf ^ 0.5i?e^/2 and k^^'^'/kf 0.8kfRe/Ha. Finally, we show that the set of least dissi- 
pative modes gives a relevant prediction for the scale of the first three-dimensional structure 
to appear in a forced, initially two-dimensional turbulent flow. This completes our numer- 
ical demonstration that the least dissipative modes can be used to simulate both two- and 
three-dimensional low-Rm MHD flows. 

1 Introduction 

Turbulence can be described as a flow where a large number of different patterns evolve in complex 
interaction with one another. The knowledge of how much energy each of them carries at a given 
time then provides a reasonably simple statistical representation of the flow. Our purpose is to 
apply this very idea to turbulence in liquid metal flows subjected to an homogeneous external 

magnetic field, by tailoring existing spectral methods to this particular prob lem. 

Although simple, these ideas express quite closely the phenomenology behind [Kolmogorov ( 1941 )'s 



famous theory of homogeneous isotropic turbulence. Here, the patterns are isotropic vortices sorted 
in three categories, according to their size Ik (or wavelength k): the large scales where energy is 
injected in the flow through some unspecified forcing, the inertial range, where mid size vortices 
pass on energy to smaller scales and the smallest scales of size fvi<^ Re^^^ where kinematic energy 
is dissipated by viscous friction {Re = UL/u stands for the Reynolds number built on velocity U 
and length L, that are typical of the large scales, as well as the fluid kinematic viscosity u). This 
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early picture has been a lot furt her refined since the n, to account for more complex effects such as 
intermittency (see Frisch ( 1995 ) or Davidson ( 2QQ4I ) for an overview). 

The description of the flow in terms of patterns is also well reflected in the more mathematical 
spectral approach of turbulence, in which the solution is sought as a decomposition over the 
elements of a basis that spans the functional space it evolves in: 



u = y^Qft)ui(x). 



(1) 



The spatial dependence (x) representing the flow patterns is carried by while the time depen- 
dence (t) appears in the coefficients of the expansion Ci only, so when ([T]) is injected into the set 
of Partial Differential Equations that gove rns the problem, th e latter reduces to a simpler sys- 
tem of Ordinary Differential Equations (see Canuto et al. ( 20061 ) for a detailed account of spectral 
methods in fluid mechanics). Apart from clear advantages in terms of simplicity and precision, 
spectral methods can also be tailored to the physical reality they describe by choosing a basis (u^) 
that represents realistic flow patterns. This basis can be obtained from the set of eigenvectors and 
adjoint eigenvectors of the operator derived from the linear part of the motion equations, with 
the boundary conditions of the problem. In incompressible homogeneous turbulence in a spatially 
periodic domain, the corresponding operator is the self-adjoint Stokes operator. Its eigenvectors 



are Fourier functions (|Constantin et al\ (|l985l )). which are classically related to vortices of wave- 
vector k. When the flow is isotropic, vortices of all shapes are present in statistically equal number, 
so they are only sorted according to their size ||k||, which facilitates the direct comparison with 
Kolmogorov's phenomenology. 



The picture is quite different for turbulence in liquid metals, where the application of a strong 
magnetic field B breaks isotropy. The fluid motion induces eddy currents that produce strong Joule 
dissipation and interact with the magnetic field to yield the Lorentz force. When the magnetic 
Reynolds number Rm is small, as in most experiments at the laboratory scale, the magnetic field 
induced in turn by these currents can be neglected so the total magnetic field is externally imposed 
and not altere d by the fluid motion. In the frame of this so-called Low Rm approximation (see 
RobertsI (1963)), the Lorentz force mainly damps velocity variations along the magnetic field lines 
so vortices tend to be elongated in this direction, resulting in a strongly anisotropic flow. This 
effect is counteracted by inertial effects that tend to break up long vortices and promote isotropy 
in the flow. Just how isotropic the flow is, is determined by the ratio between the Lorentz force 
and inertia, expressed by the interaction parameter S = aB'^L/{pU)^ where a and p are the fluid's 
electric conductivity and density. For large in a three-dimensional cubic periodic domain, when 
all vortices extend from one boundary to the other, the flow is perfectly two-dimensional, so a 
transition exists between two- an d three-dimensional turbule nce. These effects were pointed out in 
the 1960's ([MoffattJ (|l967h ) while[s ommeria fc Moreaul ( Il982l ) analysed the conditions for a channel 
flow p erpendicular to the magnetic fleld B to be quasi two-dimensional. More recently, iDavidson 
(|l997h explained how vortices evolve using the conservation of angular momentum. 
Spectral methods have been numerically implemented to study this t ype of flow in thre e-dimensional 
periodic domains in several important pieces of work, starting with Schumann ( 19761 ) who showed 
that the free decay of initially isotropic turbulence under the influence of an homogeneous mag- 
netic fleld in a three-di mensional periodic box at high S could lead to a two-dimensional state. 
IZikanov & ThessI ( 1998h found that initially isotropic MHD flows held steady on average by ap- 
plication of a forcing localised in a spherical shell of the Fourier space exhibited intermittent 
shifts between two and three-dimensional states for 5 ~ 1. Intermittency was also observed by 
[rh ess fc Zikanovl (|2007l ) in both forced and decaying MHD flows in a tri-axial ellipsoid. Most of 



these studies, however, have used the basis derived from the Stokes operator, and analysed the 
flow in terms of the modulus of the structure's wavevector /c, when clearly, anisotropy imposes 
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that vortices of same k but oriented along or across the magnetic field should undergo very dif- 
ferent levels of Joule dissipation and eventually carry very different levels of energy. Also, since 
no clear MHD equivalent to the Kolmogorov laws had been derived at the time, Kolmogorov laws 
themselves were used to impose a global cutoff frequency on k when once again, the resolution 
required to resolve the fiow completely would be expected to decrease when spanning directions 
from across to along the magnetic field direction. Therefore, determining a more " MHD-suitable" 
basis, and obtaining MHD equivalent to the Kolmogorov laws for the dissipative scales in both 
two- and three-dimensional MHD forced turbulence are the precise questions we wish to address in 
this work, by going back to the initial idea of using a basis of functions that imitates fiow patterns 
as closely as possible. We focus our attention on the configuration of a cubic domain, periodic 
in the three spatial directions, with an homogeneous magnetic field in the z direction. Although 
physically not realistic, these assumptions offer a simple but still meaningful test case for the ap- 
plication of our ideas, keeping in mind that results more directly comparable to experiments will 
have to come out of a configuration where boundaries that intercept the magnetic field lines at 
least, will be physical walls. 

In the frame of the low Rm approximation, the Lorentz force appears as a linear term in the 
Navier- Stokes equation so the linear par t of the latter is in fact the sum of the Stokes operator 
and that related to the Lo rentz force (see Roberts ( 1967h ). We have previously solved the spectral 
problem for this operator ( Potherat fc Alboussierd (|2QQ3[ )). shown that it was self- adjoint and that 
its sequence of eigenfunctions (the least dissipative modes) was able to finely mimic the anisotropic 
properties of MHD turbulence. We also showed that this sequence of modes achieved an up- 
per bound for the attractor dimension of the system that was consistent with estimates obtained 
heuristically for the size of the smallest scales. It is worth mentioning that the spectral analysis 
of the same operator, but in the case where the boundaries orthogonal to B are physical walls 
leads to a sequence of eigenfunc tions that exhibit the correct Hartn iann boundary layer profile 
in the vicinity of these walls (see Potherat fc Alboussiere (|2QQ6[ ). and iMoreaul (|199Q[ ) for a review 
of the theory of these layers). In the present work, we will therefore numerically implement our 
previously found basis in order to extract the relevant modes and determine the MHD equivalent 
of the Kolmogorov scales. In section [21 we fi rst recall and complemen t the properties of the linear 
part of the Navier-Stokes equation found in Potherat fc Alboussiere (j^03). We then implement 
this basis in an existing spectral code and determine some Kolmogorov-like laws for the small scales 
in three-dimensional MHD fiows which should serve as a criterion to resolve the fiow completely 
in section [3l Since an essential property of MHD turbulence is that it can be two-dimensional 
or three-dimensional, we devote section |4] to testing whether DNS based on the least dissipative 
modes can reproduce this feature. This leads us to find out the lengthscale of vortices in which 
three-dimensionality first appears when the intensity of the forcing is increased in an initially 
two-dimensional fiow. 



2 Principle of DNS based on the least dissipative modes 
2.1 Problem formulation 

We consider an incompressible, conducting fiuid (density p, electrical conductivity a and kinematic 
viscosity u) in a three-dimensional periodic cube Q of size Lq under imposed homogeneous and 
steady magnetic field Bgz. In the frame of the \ow-Rm approximation, the governing equations 
can be reduced to the closed system made of moment um and n iass c onser vation, which involve 
the fio w velocity u(x, t) and pressure p(x, t) only (see RobertsI (|l967l ) and Sommeria fc Moreaul 



(|l982f )). A third equation deduced from electric current conservation and the Ohm's law can be 
used to reconstruct the electric potential and the electric current a posteriori. We shall, however. 
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only need here the equations for u(x, t) and p(x, t). These can be written in non-dimensional form 
by choosing reference length L, time i^^/z^, velocity z^/I/, pressure pv'^ / L'^ and a dimensionless 
external force ||f||/L^/^ where || • || = (/ | • \^dQ.Y^'^ is the usual norm in L2{^) space. The Navier 
Stokes equations are then written: 



— u(x, t) + (u ■ V)u + Vp = V^u - Ha^V-^-^ + Gf (x, t), 
V-u = 0, 



(2) 



r3/2 



where Ha = LB. I— is the Hartmann number while G = ^^^-llfll is the Grashof number, which 

represents the forcing normalised by viscous forces (as in IPoering fc Gibbonsl ( 19951 )). Conse- 
quently, the solution of (|2]) is defined by the only two relevant control parameters Ea and 6^ in (|2]) . 
The choice of L is not straightforward as it is not imposed by the geometry. It is noteworthy that if 
it is set to L = \f^^ Xhi^v^ the governing equations depend on the single dimensionless parameter 
G jLLo? . This reference length however ignores the dynamics of the large scales present in the flow. 
One would instead expect a better suited reference length to follow the forcing scale to some extent. 
Since, however, the latter is not specified at this stage, we shall choose L = Lq, as it represents 
de facto the largest achievable scale in our problem, and denote i7ao, the Hartmann number built 
on I/Q. It is worth stressing that we shall not try to minimise or ignore the effect of the bound- 
aries where periodic conditions are applied. In particular, we shall also analyse two-dimensional 
flows where structures extend across the whole domain in the z direction. Although clearly not 
experimentally achievable, this configuration has often been used as an interestin g toy-model for 



the st u dy of the transition between two-dimensional and three-dimensional flows (jNakauchi et i 
(|l992h : [Zikanov fc ThessI (|l998[ ): iThess fc Zikanm] (|2QQ7h ). Therefore, contrarily to many previous 



studies of turbulence where periodic domains are used to represent a small volume taken out of 
an homogeneous flow, and where structures of the size of the domain should therefore be avoided, 
the conditions under which structures extend over the full domain along z will be of interest in 
this work. For this reason, the length Lq will be a meaningful parameter of the problem, wherever 
such two-dimensional vortices are considered (in section |4]). 



Two further non-dimensional numbers can be defined that are traditionally used in MHD 
turbulence: the usual Reynolds number Re = ^^^"^ , with integral length scale 

CO 

^int = ^^ / mV^E(k)dk (3) 

2 u J 

" " 

gives a measure of the intensity of turbulence (Here, k is the three-dimensional wavevector that 
appears in the Fourier transform of u, E(k) is the spectral power density of all wavevectors of 
norm k and U = {J EdkY^'^ is a reference velocity). Also, the magnetic interaction parameter 
S = cfB'^Lq/ {pU) represents the ratio of the Lorentz force to inertia. In freely decaying turbulence 
where boundaries are ignored, taking as a reference velocity from the initial velocity field and Lint 
as a reference length, S becomes the only non-dimensional parameter that governs the problems. 
In our case however, only G and Ha^ are known a priori. In this sense, they are the control 
parameters for this problem. 

The problem is fully defined by the addition of periodic boundary conditions 



u(x, t) = \x{x + a, t) 
= u(x, y -\- z^ t) 
= u{x^y^ z -\- c^t)^ a^b^cGZ 



(4) 
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and of the initial condition 

u(x,?/,z,0) = Ui{x,y,z). (5) 

These, together with the mass conservation, which simply implies that u is a solenoidal vector 
field, are taken into account by specifying that the solution u is sought in the functional space 
a solenoidal subspace of Hilbert space H^. Since the spectral method we wish to implement is 
derived from the spectral properties of governing equations, these ought to be written in abstract 
form, with help of the Helmholtz decomposition: 



Details of the mathematical framework can be found in Dymkou fc Potherat ( 2QQ9I ). The advantage 



of this form is that it gathers the linear part of the equations into a single operator that operates 
in V'^ onto itself: 

DHa=p{'^'-Ha'V-'^^ : ^ V\ (7) 

P denotes the orthogonal projection onto the subspace of solenoidal fields, and nonlinear terms 
are represented by the bilinear operator 5(u, u) = P(u • V)u. 

In the absence of magnetic field, Hao = and the system reduces to the usual Navier-Stokes 
equation. Periodic boundary condition s then ensure that the eigenfunctions of the Stokes operator 



form a basis of V'^ (jFoias et all (|2QQlh ). They can thus be used for the spectral decomposition in 
order to reduce the problem to a simpler system of ordinary differential equations. For Hao 7^ 0, the 
physical relevance of the linear part can be seen by noticing that the Lorentz force only appears 
in Dffd^- The spectral properties of this operator are therefore expected to express the mode- 
selecting dissipation that results from its action on the flow. This makes the set of eigenfunctions 
of Dffd^ a good candidate for the choice of the basis of modes requi r ed in the solution's expansion 
([1]). We have previously found these in ( Potherat fc Alboussierd ( 2QQ3[ )) and shown that they 



constituted a basis of V^, so we shall now summarise and extend these results derived from the 
spectral characteristics of the dissipation operator Djjr^. 

2.2 Spectral properties of the Dna operator for any given Ha 

Dj£^ is a linear operator. The boundary conditions are accounted for in the definition of the domain 
of the operator, defined as D{A) = V'^{Q). Since Q is bounded, the natural injectio n of V'^ into 



L2{ft) is compact, thus D^Jq^-, as an operator in ^2(1^), is compact (Foias et al\ (|2QQll )). Also, this 



operator is self- adjoint and therefore possesses a discrete set of eigenvalues (Ak) and eigenfunctions 
^ that form an orthonormal basis of the L2(yt) space. We have shown in lPotherat fc Alboussiere 



(|2QQ3[ ) that the eigenfunctions v*^ = are a subset of the usual Fourier space: 

= F,e^'2^*^-^, (8) 

with wavenumbers k = {kx^ ky^ ky) G Z^, constants G C and where j is the imaginary unit. The 
corresponding eigenvalues are 

Ak = -A^\kl + kl + kl) - Ha \ ''j (9) 
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We denote the set of all eigenvalues ([9]) by aoo{Dj^^). Since Ak represents the linear decay rate 
of mode Vk by Dna^ and Ak < 0, (Ak) and can be arranged by growing dissipation. This 
singles out Ak as a spectral parameter that naturally reflects the effects of the Lorentz force. From 
the definition (|9j), we see that for Ha = 0, |Ak|/(27r)^ reduces to the square length k'^ = ||k|p 
of the wave vector k which is the usual spectral parameter in non-MHD isotropic turbulence 
(see Figur^Ha)). In the MHD case, different values of the magnetic field B or of the reference 
length L that enter the definition of Ha yield different sets of eigenvalues (see Figur^a)-(d)). 
Such dependency is absent in the usual Fourier basis ordered by growing ||k||. The main novelty 
introduced by using this basis thus doesn't reside in the elements of the basis themselves but 
rather in the fact that they are ordered by growing values of |Ak| instead of by growing k. This 
earns these modes their denomination of least dissipative. Furthermore, we previously showed 
(jPotherat & Alboussiere (2003)) that the set of least dissipative modes required to describe the 
flow possessed the anisotropy properties predicted heuristically for such MHD flows. In the light 
of ([9]), the sequence (-Ak)^/V(27r) therefore appears as an anisotropic generalisation of the usual 
/c-sequence, and the spectral decomposition ([T]) of u can now be rewritten as 

U(x,t)= ^ CA,(t)vA,(x), (10) 

|Ak|<|A— I 

where Ck(t) are the expansion coefficients, VAk(x) are the eigenvectors of for eigenvalue Ak 

and A"^^^ defines the maximum resolution required to resolve the flow completely. 



2.3 Choice of the set of least dissipative modes 

At this point, we still lack two parameters to be able to choose the set of modes to fully resolve a 
given flow, defined by the values of G (or Re) and Hao. Firstly, the 'shape' of the set of modes is 
determined by the value of Ha only. We have however defined Hao using the domain size Lq, as a 
reference length. Clearly, for Ha to reflect the actual physics of the flow, another reference length 
L should be found that accounts for the forcing scale in one way or another. Secondly, the number 
of modes N required to fully resolve the flow or, equivalently, the largest value of |A|, |A"^^^| in 
([To]) must be determined in such a way that the flow is fully represented by its projection onto the 
set of N least dissipative modes defined by |A| < |A(A/')| = |A"^^^|. For this, the global attractor 
of the motion equations has to be entirely included in the functional subspace spanned by the 
least dissipative modes. Consequently, if dM is the dimension of this attractor, or equivalently 
the number of degrees of freedom of the flow, we must have |A™^^| > \X{dM)\- Unfortunately, 
it is difficult to obtain a precise estimate for dM- Its physical interpretation, however, can be 
easily understood: in both the non-MHD and the MHD case, the reason why dM is finite is that 
viscous dissipation introduces a cutoff at the small scale s, beyo nd which flow structures carry 



a vanishingly small amount of energy. IConstantin et al\ (|l985f ) give an elegant illustration of 



the physical meaning of these mathematical concepts. This cutoff wavelength can be estimated 
heuristically, which, in turn leads to scalings for N. The most famous example is that of the 
three-dimensional non-MHD case, where the heuris tic Kolmogorov sca le k^^^ = k^ C^Re^^^ 
(= |A"^^^|^/^/27r in our notations, and where > 1 . iKolmogorovl ( 194lh ) gives an estimate that is 



precise enough to be used as a criterion to fix the number of determining modes as A" :^ C^Re^^^ in 
a Fourier-based DNS . In tw o-dimensional turbulence, a precise estimate for the attractor dimension 
([Doering & Gibbons" (l995|)) and a heuristic scaling for the size of the smallest, or Kraichnan scales^ 



([Kraichman (1967 ): .Qhkitariil (|l989h ) coincide precisely with fc"^^^ = kk ^ G^/^(l-Mog G)^/^ where 



^max ^ |;ymax|l/2/27r . 

In the MHD case, viscous diss ipation still determine s the cutoff scale, even though Joule d issipation 
extracts energy at all scales. lAlemanv et al\ ( 1979[ ) and IPotherat fc Alboussier3 ( 2003 ) used this 
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idea, further assumed that the anisotropy k± / kz was scale-independent and that inertia balanced 
the Lorentz force at all scales to derive some heuristic scalings for the cutoff value A"^^^ and 
when S ^ 1: 

N cCoj^, (11) 



27rkf 



CxRe''. (12) 



We have here expressed A"^^^ with respect to the largest forcing scale in the problem Lf = Lo/kf 
to reflect the fact that for spatially periodic domains, the forcing scale is a relevant large scale 
that determines the small scales while the size of the computational domain isn't. Since the set of 
least dissipative modes is a subset of that of Fourier modes, these scalings can be more classically 
expressed in terms of the smallest scales across (subscript ±) and along the magnetic field by virtue 
of the properties of ([9]): 

We have b een able to partly confirm thes e scalings by finding an upper bound for the attractor 



dimension (jPotherat fc Alboussierd (|2QQ3f )). Co or Cx however remain to be evaluated, so no 
practical criterion currently exists for the number of determining modes in flows where a magnetic 
field is present. The next section is therefore devoted to searching numerically the values of L and 
^max particular, we shall estimate the lowest values of Cx for which the flow is fully resolved 
for > 1. When 6* >> 1, the flow becomes two-dimensional so the set of least dissipative modes 
becomes the two-dimensional isotropic set of Fourier modes defined hy N > k^ (or A"^^^ = (27r/cfe)^). 
When iS* << 1, the effects of the Lorentz force become small and the set of least dissipative modes 
differs little from that of the usual three-dimensional isotropic set of Fourier modes N C^Re^^"^ 
(or A"^^^ = (27rfe^)2). 

At this point, it is important to notice that a flow described by the set of least dissipative modes 
with A"^^^ determined by the rules above is resolved exactly, without any approximation, as all 
energy and dissipation containing modes are contained in the attractor. In particular, a clear 
distinction should be made between solving the equations by projection on the full set of least 
dissipative modes, which is a type of Direct Numerical Simulation, and approaches such as Large 
Eddy Simulations where part of the spectrum is modelled and not resolved. Both approaches could 
even be combined to achieve important reductions in computational cost. 

3 Determination of the exact set of modes required to re- 
solve the flow for S > 1 

3.1 Numerical system and procedure 

We base our DNS on the eigenfunctions of the dissip ation operator. Sin c e thes e are a subset of the 
usual Fourier modes, we use the code developed by iKnaepen fc MoinI ( 2QQ4I ) and IVorobev et al 



([20051) where the problem formulated in section [2Jl was implemented and fully tested. It relies on 
traditional spectral methods based on a Fourier decomposition, wit h Fast Fourier Tran sform and a 
fourth-order low- storage time-integration Runge-Kutta scheme (see Rogallol ( 198ll ) and Williamson 



19801)). The a li as error resulting froni the bilinear products is removed by phase-shifting method 



Ro£a llQl (|l98lh : [Orszag fc PattersonI (|l97lh ). which allows us to retain all of the Fourier modes but 



requires eight evaluations during each time step. We adapt this code to our needs of performing 
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(a) Ha 



(b) Ha = 630 
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(c) Ha = 3140 
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(d) Ha = 3140 



Figure 1: Iso-A curves for different values of Ha. Note that all families of curves (except those for 
Ha = 0) can be scaled down to a single family in the {k±/Ha^ kz/ Ha) 

plane where k^ = ^//c^ + ky. Values of k^^^ and kf^^ are marked on arbitrary iso-A curves to 

illustrate how they are related to A"^^^. 
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calculations using set of modes that satisfy |Ak| < |A"^^^|, simply by setting unneeded modes to 
zero when required. In all calculations presented in the whole of section [3l initial velocities are set 
to zero {u{t = 0) = 0). The flow is driven by two distinct types of constant forcing f in (|6j), that 
respectively favour two-dimensional and three-dimensional structures. The two-dimensional forcing 
is applied to Fourier modes with wavevectors kf = {kfx^kfy^kfz) G {(6, 6, 0), (7, 7, 0), (9, 9, 0)} 



f2D(x 



sm{kfx27rx) cos{kfy27ry)ex + cos{kfx27rx) sm{kfy27ry)ey 



(14) 



and tends to generate a flow with no velocity component nor velocity variations in the z-direction. 
Since the numerical algorithm would not otherwise allow the solution of the problem to be three- 
dimensional at all, we add a small constant force of amplitude e = 10~^ (relative to f) in each 
ball ||k — kf|| < 2. There are several other reasons for this choice: firstly, the forcing has to 
be a combination of the set of modes used for the expansion. In this regard, a practically z- 
independent forcing can be used to simulate both two-dimensional flows (for which the effect of 
the small three-dimensional component of the forcing falls within the numerical error) and three- 
dimensional flows. The second reason is that this type of constant weakly three-dimensional forcing 
strongly resembles that obtained in liquid metal experiment s by injecting elec t ric current thoug h 
metallic electro d es en ibedded in insulating Hartmann walls ( Sommerial ( 1986[ ). Sommerial (| 19881 ). 
Delannov et al. (|l999[ )). O ur most recent exp e riments on electrically dr iven channel flows under 



transverse magnetic fields (jKlein et 



Klein fc Potherat (|2Q1Q[ )) have indeed confirmed 
the previous theoretical prediction that in such experiments, even for high values of Ha^ inertia 
induced some slight velocity variations along t he magnetic field l ines, so that three-dimensional 
vortex instabilities such as those analysed by iThess fc Zikanovl (|2QQ7l ) do not occur in strictly 
two-dimensional, or eve n strictly quasi two-dimensional flows, but rather is some weakly three- 
dimensi onal flow (^Potherat et al\ (|2QQQ[ )). which our weakly three-dimensional forcing imitates. 
Finally, IVorobev et al\ (|2QQ5[ ) have suggested that the two or three-dimensional nature of the forcing 
had no noticeable influence on the anisotropy of intermediate and small scales. This is supported by 
the properties of the least dissipative modes, as they imply that the small scales are determined by 



G, wh ich only carries the intensity and the scale of the forcing, and Ha (jPotherat fc Alboussiere 



(|2QQ3f )). To check this point further, we have performed a series of computations in the same 
conditions as thos e described above, b ut with a three-dimensional forcing. The latter was chosen 
of the ABC type ( Mininni et al\ ( 2QQ6I )) so as to act on the three components of the velocity, and 
expressed as: 



f3D(x,t) = {cos{kfyy) + l.ls'm{kfzz))ex + {l.lcos{kfzz) ^ 0.9 s'm{kfxx))ey 
+ (0.9 cos{kfxx) + sm{kfyy))ez with kj = (6, 6, 6). 

All calculated cases are summarised in table [H 



(15) 



3.2 Determination of the length scale Lopt 

We first address the problem of choosing the best suited reference length L that enters the definition 
of the Hartmann number Ha, for a flow at given Hao and G (or Re). This problem appears only 
in three-dimensional flows as in two-dimensional flows, the least dissipative modes reduce to the 
isotropic set of two-dimensional Fourier modes. At this point, one should remember that the choice 
of the basis is arbitrary and should not have any impact on the final solution, as long as its elements 
can be combined to obtain all the energy and dissipation-carrying modes. In the particular case of 
a basis of least dissipative Fourier modes, this gives us the freedom to leave L as a free parameter 
a priori, and to fix it so as to obtain a basis that contains the least possible non-energetic, non- 
dissipative modes, that are superfluous for the description of the solution. How this can be done 
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Figure 2: Contours of spectral density of energy E{k±,kz) (colours) with iso-k and iso-A curves 
(solid lines) for several values of L at Hao = 80 and = 2.94 x 10^. 
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Figure 3: Le/t* Variations of T^Exi^) ^^e cases listed in table [TJ The minima indicate L = Lopt- 
Symbols are those from table [H Right: variations of Lopt with Hao and G for Hao = 80 and 2D 
forcing (dash-dot), Hao = 400 and 2D forcing (dashed) , Ha = 1000 and 2D forcing (dotted) and 
Ha = 400 and 3D forcing (solid). 



can be understood by analogy with the non-MHD case where the flow is expected to be isotropic 
in regions of the Fourier space located far enough from the forcing. There, the energy of a given 
mode k is expected to depend on ||k|| only. Similarly, for the spectral parameter A to be physically 
relevant to the MHD case we would expect each eigenmode of Djj^ of eigenvalue A located far 
enough from the forced modes kj to carry approximately the same amount of energy. The erratic 
nature of turbulent flows, however, makes it impossible to satisfy this condition exactly, so we shall 
instead look for the optimal value Lopt of L that minimises the functional: 



«e<T(D^^)k:A(k) = 



Ex 



(16) 



where o-{Dj£^) refers to the flnite set of eigenvalues of Dna for the numerical resolution considered, 
denotes the energy summed over all modes of eigenvalue A and ExCk) is the spectral energy 
density at point k of the iso-A surface. T^Ex (^) gives one possible overall measure of how strongly E 
varies over shells shaped according to the iso-A surfaces in the Fourier space. In practice, we start 
from a "traditional" DNS resolved up to the Kolmogorov scale (these cases are gathered in table [1]), 
and therefore over-resolved in the MHD case (on the basis that the attractor dimension decreases 
monotonically when Ha increases (jPotherat fc Alboussier j (|2003[ )). This yields a reference solution 



from which E{k±^kz) can be extracted. We then calculate the minimum of functional T^ExiE) 
numerically (the variations of T^Ex{E) are shown on flgure[3l left). This is illustrated on a typical 
example for Hao = 80 and G = 2.94 x 10^ on flgure [2] where the sets of iso-A curves are plotted 
for several values of L along with the contours of E{k±^kz). One sees that the iso-A curves 
corresponding to Lopt/^o = 0.3 x 27r on Figure [2](c) follow the energy distribution well, as opposed 
to iso-k lines, shown on Figure [2](a) that cross many different levels of energy. This shows that 
the basis of the least dissipative modes does carry the morphology of the energy distribution quite 
realistically, provided we choose L :^ i^opt- It can be seen from the variations of Lopt with G for 
Hao ^ {80,400} on flgure [3l that it depends little on either G or Hao^ around 0.3Lo x 27r for 
Hao = 80 and 0.2Lo x 27r for Hao = 400. The fact that it still varies a little with Hao is certainly 
due in part to the "non-universality" introduced by the forcing, as the energy distribution clearly 
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Hao 


G 


Tlx X Tly X Tlz 


Re 


Cx 




a£;(0.5) 


Haopt/'27i 


'S'opt 


symbol 


2D forcing 


80 


2.67x10^ 


064 X 064 X 064 


70 


0.55 


1.33 


0.999 


23.70 


8.02 


A 


80 


7.34x10'^ 


128 X 128 X 128 


97 


0.81 


2.07 


0.996 


27.80 


7.97 


o 


80 


1.47x10^ 


128 X 128 X 128 


159 


0.63 


1.43 


0.998 


24.80 


3.87 


V 


80 


2.94x10^ 


128 X 128 X 128 


194 


0.57 


1.23 


0.996 


24.90 


3.20 


< 


80 


3.34x10'^ 


128 X 128 X 128 


195 


0.57 


1.23 


0.995 


24.50 


3.08 


□ 


80 


6.67x10^ 


128 X 128 X 128 


264 


0.50 


0.98 


0.998 


24.20 


2.21 


o 


400 


6.67x10^ 


128 X 128 X 128 


216 


0.71 


1.13 


0.995 


65.80 


20.04 


▲ 


400 


1.00x10^ 


256 X 256 X 256 


245 


1.11 


2.07 


0.986 


77.40 


24.45 


• 


400 


1.33x10^ 


128 X 128 X 128 


282 


0.53 


0.93 


0.998 


42.70 


6.47 


T 


400 


2.00x10^ 


256 X 256 X 256 


343 


0.90 


1.61 


0.994 


61.40 


10.99 


< 


400 


6.67x10^ 


512 X 512 X 256 


575 


1.28 


1.09 


0.986 


48.25 


4.05 


■ 


1000 


1.33x10^ 


512 X 512 X 256 


935 


1.08 


1.52 


0.998 


112.70 


13.58 


A 


1000 


2.67x10^ 


512 X 512 X 512 


1140 


0.94 


1.3 


0.996 


86.00 


6.44 


• 


3D forcing 


400 


1.5x108 


256 X 256 X 256 


465 


0.71 


1.28 


1.000 


95.00 


19.41 


■ 


400 


6.0x10^ 


256 X 256 X 256 


512 


0.62 


1.19 


0.999 


71.00 


9.84 


T 



Table 1: Summary of all cases calculated with initial condition u(t = 0) = 0: Grashof number 
G, embedding spectral resolution Ux x riy x riz, Reynolds number Re, resolution {Cx and C/^), 
fraction of the total energy asn contained in modes with lower |A| than the value given by ([TTj), 
"optimal" Hartmann number Haopt = Hao{Lopt/ Lq), and "optimal" interaction parameter 5'opt = 



departs from the iso-A lines in the vicinity of the forced modes. At a given Re, or G, the influence 
of these modes increases with Hao, as for higher Hao, the energy tends to stay closer to the {k^^ ky) 
plane, which brings the smallest scales closer to the forced modes kj. For the purpose of performing 
DNS based on the least dissipative modes, a precise determination of Lopt is however not necessary 
as energy and dissipation spectra E{X) and D{X) obtained with L departing by around ±30% from 
Lopt, using only modes in the region |A| < |A"^^^| (where A"^^^ was flxed according to (p!6|) derived 
in the next section) yielded no signiflcant discrepancy with those obtained from calculations based 
on Lopt exactly. This robustness also confirms that as long as the iso-A curves follow the contours 
of energy well enough in the vicinity of the small scales, then the set of Fourier modes determined 
by A"^^^ contains very few non-relevant modes. Also, since Hao^t = BLo^^^J cr / {pu) gives the most 
physically relevant measure of the Lorentz force, we shall now prefer it to Hao to express the laws 
for the small scales (fT2]) and (p!3|) . 

3.3 Scaling laws for A"^^"" 

Having chosen L = I/opt ? we have fixed a family of modes, indexed by the corresponding sequence 
of values of A. We now need to know how many of these modes are required to resolve the fiow 
fully, for given values of Ha^^t and G (or Re). A usable estimate for this number is obtained 
through a value for the numerical constant Cx that appears in the scaling law for the smallest 
scales A"^^^(iJaopt, Re) (p^ . To find it, we select four cases covering different values of Ha, G, two 
and three-dimensional forcing. In each case, we first calculate the established state with resolution 
up to the Kolmogorov scale /c^ (summarised in table [1]). Since this case is over-resolved, it serves 
as a reference for the energy and dissipation distribution in the Fourier space. We then recalculate 
several times the same fiow, but resolved up to |A|^/^ = |A^^*|^/^ = CxRe^^'^ with different values 
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of Cx and compare the corresponding power and dissipation density spectra E{X) and D{X) to 
those obtained in the reference DNS resolved to the Kolmogorov scale. Finally, the impact of the 
reduction in resolution on time dependent-flows is assessed by applying the same procedure to 
the freely decaying flow that follows a shutdown of the forcing in the established regime in the 
reference case, at t = tdecay 

Figure H] summarises all calculated cases along with resolution and embedding resolution. The 
latter is of no incidence on the solution but gives a measure of the reduction in computational cost 
incurred by using our "A-based" approach, and this, even though the spectral code we are using 
hasn't been optimised for it. 



The time-averaged energy distributions in the {k±^ kz)-p\8ine {k± = y /c^ + ky) show no visible 

discrepancy between the reference case and those for Cx G [0.29,0.59]. This indicates that even 
with the lowest resolution, which uses up to 64 times less modes than the reference case, the 
energy distribution, and the flow anisotropy are still qualitatively well rendered. An inspection of 
the corresponding A-based energy and dissipation spectra from figure [5] and [6] confirms and refines 
the picture: the small energy and dissipation pile-up that inevitably occurs at the high-A end of 
the spectrum certainly remains confined there for all Ha for Ca ^ 0.5. It does, however tend to 
slightly spread towards the higher end of the spectrum for lower values of Ca, particularly in the 
dissipation spectra and in the cases at lower Re. Even though it is only pronounced in the Hao = 80 
case, this propagation of error toward larger scales is a usual symptom of under-resolution, and 
can be more easily spotted on the dissipation spectra. This error on the dissipation is further 
revealed when the flow is freely decaying. For each of our four reference cases, we have calculated 
such flows starting from an initial state in the established regime resolved up to the Kolmogorov 
scale. In each case, the subsequent evolution of the flow without forcing was calculated several 
times from this same initial condition, for the same maximum resolutions as those used to calculate 
the established flows. The evolution was calculated over 20 Joules times, after which the flow had 
lost most of its energy. As for the dissipation spectra in the established state, it turns out that a 
discrepancy between reference case and cases resolved with Ca < 0.5 is visible in the evolution of 
both the total energy and of the energy in the field direction. Cases resolved with Cx ^ 0.5, on 
the contrary, match the reference case to a great precision, both when the flow is estabhshed and 
freely decaying. As a matter of fact, the decay curves for Cx ^ 0.5 cannot be distinguished from 
those of the reference case on the graph. 



To quantify the precision reached for a given value of Cx = y |A"^^^|/(47r^/cji?e) over a wider range 
of parameters than those of the 4 reference cases calculated above, we define a reduced spectral 
parameter normalised by scaling ([12]): / = ^J\X\/ {An'^k'jRe), such that for / = Ca, A = A"^^"". We 

have calculated the variations of total energy T;e{1) and dissipation contained in the spectral 

subspace enclosed in the iso-A curve for each value of / < Ca for a selection of cases resolved beyond 
Cx = 0.5 (summarised in tabled] along with their resolution expressed in terms of Cx and C^). 
The results are illustrated on figure [T] Firstly, it turns out that for a given value of /, the ratio 
aE{l) of TiE{l) to the total energy T^EiCx) is constant for all calculated cases, regardless of the 
values of iJa, G and of the nature of the forcing (with, in particular, aE(l = 0.5) 0.99 no matter 
how high Cx is, as shown in table [1]). In other words, the precision attained for a given value of Cx 
remains essentially constant when Ha and G are varied beyond their values in the four reference 
cases calculated above. This brings further confirmation of the validity of scaling laws (p!2]) . and of 
their independence of the nature of the forcing. Secondly, the variations of Y^EiJ) and T,d(}) also 
comfort us in the choice of Ca — 0.5 as the minimum cutoff scale for full resolutions: this value is 
indeed located at the beginning of a plateau where further increase of resolution hardly brings any 
variation in the total energy and dissipation of the solution. Smaller values of Ca, on the other 
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Figure 4: Logarithmic energy distribution in the (/c^, /c^) -plane. Blue dots correspond to low 
energy modes. Each column represents flows calculated with the same control parameters and the 
same forcing (indicated at the top), but with different resolutions, determined by the value of Cx 
or, equivalent ly, by the spectral domain of resolution defined by A < A^^*, both indicated below 
each graph. Calculations from the first line are resolved up to the Kolmogorov scale k^ = C^Re^^^, 
with Re = 97, 245, 1140, 512 (see table [T]). The dashed lines indicate the embedding resolutions 
used in our code (in brackets). Modes in the white area within this rectangular domain are set to 
0. 
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Figure 5: Top: energy density spectrum in A-shells (statistically steady flow), middle-top: dissipa- 
tion density spectrum in A-shells (statistically steady), middle-bottom: energy density spectrum in 
/c-shells (statistically steady), bottom: evolution of the total kinetic energy of freely decaying flows 
normalised by the total energy at the time when the forcing was shut down tdecay For given Hao 
and G, initial conditions are taken from the same statistically steady reference flow resolved up 
to the Kolmogorov scale for all values of A^^*. Each column presents data from the corresponding 
cases from figure HI from which different resolutions are represented by the following curves: figure 
[3 left, right and figure\^ left: Cx = 0.35 (dash-dot) , Cx = 0.47 (dash), Cx = 0.59 (solid), figure\^ 
right: Cx = 0.29 (dash-dot) , Cx = 0.38 (dash), Cx = 0.48 (solid). Dotted lines correspond to the 
reference case resolved up to the Kolmogorov scales Cf^Re^^^ . 
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/2D, Hao = 1000, G = 2.7 X 10^ 
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Figure 6: (figure [5] continued) 
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Figure 7: Energy T^e{1) {top) and dissipation T^d{1) {bottom) contained in the subspace {A, |A| < 
Air'^PRe} vs. I. The symbols are those from table 1 and placed at / = 0.5. Left: two-dimensional 
forcing, right: three-dimensional forcing. 



hand, may fall outside this region and calculations at the corresponding resolution may thus fail to 
capture noticeable fractions of the total energy and dissipation. On these grounds we shall finally 
propose the following scaling for A"^^^: 



27rkf 



0.5Re 



1/2 



(17) 



The values of k^ ^^ and k^^^ can be directly de duced from that of A"^^^ through ([T3|) to quantify 
the scalings from IPotherat & Alboussierd (|2QQ3[ ) as: 



OMf 



Re 
Ha. 



opt 



^max 



(18) 



Also, the values of A"'"''' 
Re. Denoting Gf 



can be expressed as a function of G, which is known a priori, unlike 



where || • ||2/ represents the L2 norm for a domain of volume 



= {L/kf)^, the corresponding graph, on figure [HI suggests the scaling: 



(19) 
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Figure 8: Scaling law for the small scales expressed as a function of G. The symbols correspond 
to cases from table [TJ 



Finally, beyond the identification of A"^^^, A-based spectra from figure [5] and [6] exhibit an interesting 
feature, as a remarkable steep tail is present in both energy and dissipation A-spectra at high values 

1 /2 

of A. In all cases, it starts when A reaches the value of the eigenvalue Aj = —HaJ^^ of the first 
mode with a wavevector orthogonal to B. Since 5 > 1 in all our calculations, such modes, of 
the form ((0,A:a)), are located outside the Joule cone and therefore strongly suppressed by Joule 
dissipation. This explains why they carry very little energy. 

3.4 Practical use of the scaling laws 

The results of the present section now allow us to put forward a simple procedure to resolve three- 
dimensional MHD fiows in periodic domains: firstly, Lopt and Lint can be calculated at every time 
step as the numerical simulation progresses (as is already usual for Lint). When 5 > 1, (fT7|) or 
([T8|) then provide criteria for the resolution necessary to represent a three-dimensional MHD fiow 
completely. 

When 5* < 1, the fiow becomes progressively more isotropic and so does the set of least dissipative 
modes. Accordingly, the resolution required to fully resolve the fiow becomes higher than that 
predicted by scalings (fT7|) or (p!8|) . In this case. Lint and Lopt are still determined "on the fiy" but 
the usual Kolmogorov criterion must be used instead of (pT|) or (p!8|) . 

When 6* >> 1, the fiow can be either two or three-dimensional, which poses an important question 
about the ability of the least dissipative modes to represent the fiow accurately: on the one hand, 
when A"^^^ exceeds a value that depends on Ha only, a first three-dimensional mode appears in 
the set of least dissipative modes, independently of the behaviour of the fiow itself. When G is 
increased from 0, on the other hand, a first three-dimensional physical mode appears in the fiow 
at the actual transition between two- and three-dimensionality, independently of the method used 
to calculate it. We shall examine in the next section whether both coincide. This will tell us 
whether the least dissipative modes can be used for the simulation of MHD turbulence, regardless 
of whether it is two- or three-dimensional. 
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4 Least dissipative modes at the transition between two- 
dimensional and three-dimensional turbulence 

4.1 Two- vs. three-dimensional sets of least dissipative modes 

We now focus on the question of how to calculate flows using the least dissipative modes at the 
transition between two- and three-dimensional MHD turbulence. The set of least dissipative modes 
can either contain only two-dimensional modes or both two or three-dimensional modes, depending 
on the value of A"^^^. A transition between these two types of sets therefore occurs at a value 
j^max ^ ^^3B which the curve A = A"^^^ encloses at least one mode with kz > 1 (bold dashed 
line in Figure [T]^d)). According to our previous work ( Potherat fc Alboussierd (|2QQ3[ )) and in the 



present notations, the first three-dimensional mode in this sense is associated to the eigenvalue 



|^3D|_2^^, (20) 
' ' 27r ' 

and the modulus of the corresponding wavevector in the plane across the magnetic field lines is: 



|kf I = - 1- (21) 

It is important to notice that although a flow represented by a set comprising three-dimensional 
modes is potentially three-dimensional, it isn't necessarily three-dimensional. Instead it can be 
either two-dimensional or in a state of intermittency between the two states, as in lZikanov fc ThessI 



(| 19981 ). if the coefficients of the three-dimensional modes in expansion (fTQ|) are or intermittently 
become 0. This behaviour is determined by the flow dynamics, independently of the basis chosen 
to represent it (provided the flow is correctly resolved, obviously.). We shall now compare the first 
least dissipative three-dimensional mode to the first three-dimensional mode that appears in the 
flow. 



4.2 Numerical procedure 

We use the same numerical solver as that described in section 13.11 and also the same type of two- 
dimensional forcing (HH)- On the top of previous calculations initialised with the fluid at rest, 
we now perform two additional series of calculations, at Hao = 80 and Ha^ = 400 respectively, as 
follows: we start with fixed Hao^ low G and the fluid initially at rest. We look for a statistically 
steady two-dimensional solution and let it reach a well developed, turbulent state (after a time of 
the order of 100 - 2006'o, or dimensionally, 100-200 Joule times tj = p/{aB^)). With this latter 
state as the initial condition, we perform the next calculation by increasing the Grashof number 
by 15%, and repeat the procedure until three-dimensionality appears. 

In all simulations the numerical resolutions Ux x Uy x Uz are chosen as the smallest powers of 
2 such that the resolution domain encloses the A = 1.5A"^^^ curve and satisfies k^^^ > 1.2/cx = 
1.2G^^/^(1 + log G)^^^. This way, the flow is well resolved whether in a state of two-dimensional 
turbulence or in a state of three-dimensional MHD turbulence. Since Lopt cannot be determined in 
two-dimensional flows but varies little for a given Ha, we take the approximate values Lopt{Hao = 
80) = 24 and Lopt{Hao = 400) = 55, (see figure [3]). 

4.3 Measure of three-dimensionality 

In order to track three-dimensionality near the transition, we define two quantities to characterise 
it. The first one expresses how physical quantities depend on so we shall call it morphological 
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three-dimensionality and define it as 

\ 1/2 

«3D= I / {f{z)-lfdz\ , (22) 

where f{z) expresses the ratio between the two and three-dimensional parts of the RMS velocity 
fluctuations in the plane z = const: 

Lo J {<{M'{x,y,z')f>tf/^dxdy 

^^'^ ^ /(<(u^(x,y,z))2>,)i/2dxd^dz • ^^^^ 
n 

Here, < • >t denotes averaging with respect to time and = u— < u >t is the local velocity 
fluctuation. asD gives a global measure of morphological three-dimensionality as it expresses an 
average ratio of the three-dimensional to the two-dimensional part of the velocity fluctuations. 

The second type of three-dimensionality is expressed as the ratio of the energy in the z direction 
to that in the x and y direction. We shall therefore call it kinematic three-dimensionality: 



In theory, there is no reason for the first appearance (in the sense of growing G) of these types of 
three-dimensionality not to take pi 
name k^^^ and k^^^ respectively. 



three-dimensionality not to take place in vortices of distinct wavelength, which we shall therefore 



4.4 First three-dimensional modes and relevance of the least dissipative 
modes to transitional flows 

On the cases initialised with the fluid at rest, we flnd that both asD and f3sD jump to flnite 
values at the same value of the forcing G^^{Haopt). By contrast, when the forcing is increased 
progressively, morphological three-dimensionality appears at a lower critical value of G than dy- 
namical three-dimensionality. We have identifled /c^^" and k^^ by calculating the quantities 
E^'^{k±) = E±{k±,kz) and E^^{k±) = ^z{k±,kz) respectively. Both are plotted on 

flgure [9] for the flrst value of the forcing where three-dimensionality was observed. These quan- 
tities indeed remain at noise level for two-dimensional flows. When morphological {resp. kine- 
matic) three-dimensionality appears, several peaks rise in the proflle Ef''{k±) {resp. E^^{k±)) 
at k± = /c^^" {resp. k± = k^^^). Further peaks also appear around /c^^" and k^^^ . This is 
due to the fact that three-dimensionality can only be detected in slightly supercritical regime. 
Furthermore, since the maximum of the iso-A curve in (/c^, kz) is not only very "flat" but can also 
be located at a non-integer value of /c^, several peaks are expected to rise around the maximum. 
This is all the more true at high Ha. Keeping this in mind, one still sees that at the lowest forc- 
ings where either morphological or kinematic three-dimensionality were detected, both appeared in 
columnar vortices of approximately the same wavelength k^^^ k^^^ . Importantly, this value is 
consistent with the theoretical estimate (|2T]) for /c^^, albeit a little smaller in the case Ha^ = 400. 
On the top of the iso-A curve being very flat at Ha^ = 400, this shift towards larger scales can be 
explained by the fact that the higher iJa, the higher the value of G at which three-dimensionality 
appears, and the higher the turbulence intensity when this happens. In two-dimensional turbu- 
lence, inertial transfer increases the energy of the large scales, that are therefore more prone to 
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Figure 9: Profiles of Ef^{ki_) and Ef^{k±) for Hao = 80 (left) and Hao = 400 (right): the 
corresponding flows are weakly three-dimensional. Curves marked with 'o' symbols indicate cases 
initialised in a stabilised state at slightly lower forcing while curves without them correspond to 
flows initialised at rest. The vertical dashed lines mark the theoretical values of given by (|2T]) . 



exhibit instabilities leading to the appearance of three-dimensionality. Among the least dissipative 
modes that dissipate energy at about the same rate, this favours those with the larger scales, over 
the strictly least dissipative one predicted by (|2T]) . 

Importantly, one sees on figure [9] that k^^^ o:^ k^^^ :^ k^^ is independent of the flow's initial 
conditions, even though a^^D and P^^d aren't. In other words, even in cases where morphological 
and dynamical three-dimensionality appear successively (in the sense of growing G) they do so 
in vortices of the same transverse wavelength (|2T]) . This implies that one can use the set of least 
dissipative modes together with scalings (pT|) or ([18]) in order to determine a priori the exact set of 
modes required to resolve both transitional and three-dimensional flows completely, provided Lopt 
is known (It can be obtained from the calculation of a three-dimensional flow at the same value of 
i^ao, for instance.). For flows that lay at the transition between two- and three-dimensionality, a 
slight over-resolution is advisable that will absorb the peaks of three-dimensionality that appear 
around k^^ ^ kl^^ k^^^ . 

It is quite remarkable that for the forcing (and the forcing scale) we have chosen, k^^ follows 
(|2T]) rather well. Just how universal this behaviour is, however, remains to be clarified. For a 
sufficiently turbulent two-dimensional flow forced at /c/ > /c^^, the inverse energy cascade can 
be expected to transfer energy back to k^^ where three-dimensional vortices would form. More 
generally, our recent experiments on MHD turbulence in cubic box have shown that the appearance 
of three-dimensionality was go verned by a subtle i nterp lay between inertia and the Lorentz force 
at the scale of each structure ( Klein fc Potherat ( 2010[ )). The former is determined on the one 



hand by the forcing, which arbitrarily injects energy in the flow and, on the other hand, by the 
turbulent redistribution of energy amongst structures. Flows where turbulence is absent or too 
weak to sufficiently erase the non-universal trace of the forcing, therefore don't exhibit the ideal 
behaviour predicted by (|2T]) . This was spectacularly illustrated in our experiment where at low Ha 
and low i?e, the destabilisation of a periodic array of columnar vortices led to remarkable steady 
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three-dimensional Y-shaped vortices. 



5 Conclusions 

In this article, we have shown that DNS of Low-Rm MHD turbulence in a three-dimensional pe- 
riodic domain could be achieved by using the sequence of least dissipative eigenmodes from the 
dissipation operator instead of the traditional Fourier basis. Not only is this technique far more 
cost effective at fully resolving the flow without modelling, but it also enlightens some of its prop- 
erties that don't appear otherwise. Indeed, the iso-energy lines follow the lines of constant linear 
decay rate A well in regions of the spectral space that are not directly influenced by the forcing. 
Furthermore, energy and dissipation spectra expressed in terms of the eigenvalue A associated to 
these modes instead of /c, exhibit a clear cutoff that identifies modes located inside the Joule cone, 
and therefore strongly suppressed by Joule dissipation. Most importantly, analysing this spectra 
for 5 > 1 allowed us to derive laws that play the role of Kolmogorov laws, of determining the small 
scales in MHD turbulence: V^A^/(27rfc/) 0.5i?e^/^ or v1^^/(27r/c/) ^ 0A7G^-'^^. Finally, 
MHD flows in a periodic domain can be resolved as follows: Lopt and Lint can be obtained on 
the fly, by minimising functional T^Ex every time step (see section 1X2]) . The discrete sequence 
of values of A then follows from (|9]), and ultimately, the small scales are obtained using our new 
scalings (fT7|) if ^opt ^ 1, or the Kolmogorov laws if 6* < 1. 

In the last part of this work, we also showed that the set of least dissipative modes encompassed 
the modes that first exhibit three-dimensionality when the forcing was increased from either zero 
or from that of a two-dimensional flow. This proves that the set of least dissipative modes is also 
suitable for the resolution of transitional flows, and not only for three-dimensional flows. On the 
top of this, for two-dimensional flows, that occur in the limit of large the Lorentz force vanishes 
so the set of least dissipative modes coincides with the usual set of two-dimensional Fourier modes. 
They can therefore be used in conjunction with Kraichnan's law for the size of the smallest scales 
|;ymax|i/2/(27r) - G^^^. The Icast dissipative modes can therefore be used to calculate MHD flows 
in a periodic box for all values of S. 



Finally, we wish to underline the large potential field of application of the method presented in 
this work. The initial idea was to use a basis of modes that already incorporates the main consti- 
tutive structures of the flow, so as to save the costs of having to reconstruct them using elements 
of a less suited basis. In the present case, the basis of least dissipative modes readily rendered the 
anisotropic properties of MHD turbulence. Using this basis therefore reduced the cost of DNS by 
confining the spectral domain of resolution to that strictly relevant to the flow dynamics. This 
procedure can clearly be extended to MHD and non-MHD problems with more complex boundary 
conditions. We have recently shown that the orthogonal set of least dissipative modes in a channel 
flow with transverse magnetic field were exponential functions that incorporated the profi l e of t he 
very thin Hartmann boundary layers which arise along the walls ( Dymkou fc Potherat ( 2QQ9h ). 
Currently, channel fiow DNS are limited to Ha below a few hundred because of the computational 
cost involved in meshing these layers. Using the least dissipative modes for this problem not only 
brings the same benefits as in the periodic case studied in the present work, but it also eliminates 
the difficulty posed by the Hartmann layers as they do not have to be reconstructed nor meshed. 
As a spectacular consequence, the computational cost of DNS based on these modes decreases 
with Ha instead of increasing as in current methods based on Tchebychev Polynomials. Using the 
least dissipative modes is therefore not only beneficial to the simulation of turbulent fiows but also 
potentially to all fiows where the reconstruction of anisotropic structures with unsuited elements 
incurs computational costs far beyond those strictly required by the dynamics. 
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